﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace OseaCS
{


    public static class QrsDetH
    {
        static public int SAMPLE_RATE = 200;	// Sample rate in Hz.
        static public double MS_PER_SAMPLE = ((double)1000 / (double)SAMPLE_RATE);
        static public int MS10 = ((int)(10 / MS_PER_SAMPLE + 0.5));
        static public int MS25 = ((int)(25 / MS_PER_SAMPLE + 0.5));
        static public int MS30 = ((int)(30 / MS_PER_SAMPLE + 0.5));
        static public int MS80 = ((int)(80 / MS_PER_SAMPLE + 0.5));
        static public int MS95 = ((int)(95 / MS_PER_SAMPLE + 0.5));
        static public int MS100 = ((int)(100 / MS_PER_SAMPLE + 0.5));
        static public int MS125 = ((int)(125 / MS_PER_SAMPLE + 0.5));
        static public int MS150 = ((int)(150 / MS_PER_SAMPLE + 0.5));
        static public int MS160 = ((int)(160 / MS_PER_SAMPLE + 0.5));
        static public int MS175 = ((int)(175 / MS_PER_SAMPLE + 0.5));
        static public int MS195 = ((int)(195 / MS_PER_SAMPLE + 0.5));
        static public int MS200 = ((int)(200 / MS_PER_SAMPLE + 0.5));
        static public int MS220 = ((int)(220 / MS_PER_SAMPLE + 0.5));
        static public int MS250 = ((int)(250 / MS_PER_SAMPLE + 0.5));
        static public int MS300 = ((int)(300 / MS_PER_SAMPLE + 0.5));
        static public int MS360 = ((int)(360 / MS_PER_SAMPLE + 0.5));
        static public int MS450 = ((int)(450 / MS_PER_SAMPLE + 0.5));
        static public int MS1000 = SAMPLE_RATE;
        static public int MS1500 = ((int)(1500 / MS_PER_SAMPLE));
        static public int DERIV_LENGTH = MS10;
        static public int LPBUFFER_LGTH = ((int)(2 * MS25));
        static public int HPBUFFER_LGTH = MS125;

        static public int WINDOW_WIDTH = MS80; // Moving window integration width.
        static public int FILTER_DELAY = (int)(((double)DERIV_LENGTH / 2) + ((double)LPBUFFER_LGTH / 2 - 1) +
                                               (((double)HPBUFFER_LGTH - 1) / 2) + PRE_BLANK);  // filter delays plus 200 ms blanking delay
        static public int DER_DELAY = WINDOW_WIDTH + FILTER_DELAY + MS100;

        //TODO: pre blank nao pertence mais ao escopo da QRSDET e sim da QRSDETH
        static public int PRE_BLANK = MS200;
    }

    
    unsafe public static class QrsDet
    {

        /*****************************************************************************
        FILE:  qrsdet.cpp
        AUTHOR:	Patrick S. Hamilton
        REVISED:	12/04/2000
          ___________________________________________________________________________

        qrsdet.cpp: A QRS detector.
        Copywrite (C) 2000 Patrick S. Hamilton

        This file is free software; you can redistribute it and/or modify it under
        the terms of the GNU Library General Public License as published by the Free
        Software Foundation; either version 2 of the License, or (at your option) any
        later version.

        This software is distributed in the hope that it will be useful, but WITHOUT ANY
        WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
        PARTICULAR PURPOSE.  See the GNU Library General Public License for more
        details.

        You should have received a copy of the GNU Library General Public License along
        with this library; if not, write to the Free Software Foundation, Inc., 59
        Temple Place - Suite 330, Boston, MA 02111-1307, USA.

        You may contact the author by e-mail (pat@eplimited.edu) or postal mail
        (Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville,
        MA 02143 USA).  For updates to this software, please visit our website
        (http://www.eplimited.com).
          __________________________________________________________________________

        This file contains functions for detecting QRS complexes in an ECG.  The
        QRS detector requires filter functions in qrsfilt.cpp and parameter
        definitions in qrsdet.h.  QRSDet is the only function that needs to be
        visable outside of these files.

        Syntax:
	        int QRSDet(int ecgSample, int init) ;

        Description:
	        QRSDet() implements a modified version of the QRS detection
	        algorithm described in:

	        Hamilton, Tompkins, W. J., "Quantitative investigation of QRS
	        detection rules using the MIT/BIH arrhythmia database",
	        IEEE Trans. Biomed. Eng., BME-33, pp. 1158-1165, 1987.

	        Consecutive ECG samples are passed to QRSDet.  QRSDet was
	        designed for a 200 Hz sample rate.  QRSDet contains a number
	        of static variables that it uses to adapt to different ECG
	        signals.  These variables can be reset by passing any value
	        not equal to 0 in init.

	        Note: QRSDet() requires filters in QRSFilt.cpp

        Returns:
	        When a QRS complex is detected QRSDet returns the detection delay.

        ****************************************************************/

        //#include <mem.h>		/* For memmov. */
        //#include <math.h>
        //#include "qrsdet.h"
        //#define PRE_BLANK	MS200

        //// External Prototypes.

        //int QRSFilter(int datum, int init) ;
        //int deriv1( int x0, int init ) ;

        //// Local Prototypes.

        //int Peak( int datum, int init ) ;
        //int median(int *array, int datnum) ;
        //int thresh(int qmedian, int nmedian) ;
        //int BLSCheck(int *dBuf,int dbPtr,int *maxder) ;

        //int earlyThresh(int qmedian, int nmedian) ;

        public static double TH = 0.475;
        public static int[] DDBuffer = new int[QrsDetH.DER_DELAY]; /* Buffer holding derivative data. */
        public static int DDPtr;	
        public static int Dly = 0;
        public static int MEMMOVELEN = 7 * sizeof(int);


        public static int det_thresh, qpkcnt = 0;
        public static int[] qrsbuf = new int[8], noise = new int[8], rrbuf = new int[8];
        public static int[] rsetBuff = new int[8];
        public static int rsetCount = 0;
        public static int nmedian, qmedian, rrmedian;
        public static int count, sbpeak = 0, sbloc, sbcount = QrsDetH.MS1500;
        public static int maxder, lastmax;
        public static int initBlank, initMax;
        public static int preBlankCnt, tempPeak;


        unsafe public static int QRSDet(int datum, int init)
	    {
            int fdatum, QrsDelay = 0;
            int i, newPeak, aPeak;

            //	Initialize all buffers to 0 on the first call.	

            if (init != 0)
            {
                for (i = 0; i < 8; ++i)
                {
                    noise[i] = 0; // Initialize noise buffer 
                    rrbuf[i] = QrsDetH.MS1000; // and R-to-R interval buffer. 
                }

                qpkcnt = maxder = lastmax = count = sbpeak = 0;
                initBlank = initMax = preBlankCnt = DDPtr = 0;
                sbcount = QrsDetH.MS1500;
                QrsFilt.QRSFilter(0, 1); // initialize filters. 
                Peak(0,1);
            }

            fdatum = QrsFilt.QRSFilter(datum, 0); // Filter data. 


            // Wait until normal detector is ready before calling early detections. 

            aPeak = Peak(fdatum, 0);

            // Hold any peak that is detected for 200 ms
            // in case a bigger one comes along.  There
            // can only be one QRS complex in any 200 ms window.

            newPeak = 0;
            if (aPeak != 0 && preBlankCnt == 0)		    // If there has been no peak for 200 ms
            {								            // save this one and start counting.
                tempPeak = aPeak;
                preBlankCnt = QrsDetH.PRE_BLANK;	    // MS200
            }
            else if (aPeak == 0 && preBlankCnt != 0)    // If we have held onto a peak for
            {									        // 200 ms pass it on for evaluation.
                if (--preBlankCnt == 0)
                    newPeak = tempPeak;
            }
            else if (aPeak != 0)					    // If we were holding a peak, but
            {									        // this ones bigger, save it and
                if (aPeak > tempPeak)		    	    // start counting to 200 ms again.
                {
                    tempPeak = aPeak;
                    preBlankCnt = QrsDetH.PRE_BLANK;    // MS200
                }
                else if (--preBlankCnt == 0)
                    newPeak = tempPeak;
            }

            	
            //newPeak = 0 ;
            //if((aPeak != 0) && (preBlankCnt == 0))
            //    newPeak = aPeak ;
            //else if(preBlankCnt != 0) --preBlankCnt ;
            

            // Save derivative of raw signal for T-wave and baseline shift discrimination. 
        	
            DDBuffer[DDPtr] = QrsFilt.deriv1(datum, 0);
            if (++DDPtr == QrsDetH.DER_DELAY)
                DDPtr = 0;

            // Initialize the qrs peak buffer with the first eight 	
            // local maximum peaks detected.						

            if (qpkcnt < 8)
            {
                ++count ;
                if(newPeak > 0) count = QrsDetH.WINDOW_WIDTH;
                if(++initBlank == QrsDetH.MS1000)
                {
                    initBlank = 0;
                    qrsbuf[qpkcnt] = initMax;
                    initMax = 0;
                    ++qpkcnt;
                    if (qpkcnt == 8)
                    {
                        qmedian = median(qrsbuf, 8);
                        nmedian = 0;
                        rrmedian = QrsDetH.MS1000;
                        sbcount = QrsDetH.MS1500 + QrsDetH.MS150;
                        det_thresh = thresh(qmedian, nmedian);
                    }
                }
                if (newPeak > initMax)
                    initMax = newPeak;
            }
            else	// Else test for a qrs. 
            {
                ++count;
                if (newPeak > 0)
                {
        			
        			
                    /* Check for maximum derivative and matching minima and maxima
                       for T-wave and baseline shift rejection.  Only consider this
                       peak if it doesn't seem to be a base line shift. */
        			   
                    if (BLSCheck(DDBuffer, DDPtr, ref maxder) == 0)
                    {
                        // Classify the beat as a QRS complex
                        // if the peak is larger than the detection threshold.
                        if (newPeak > det_thresh)
                        {
                            itaShifRight(rrbuf, MEMMOVELEN); //memmove(&qrsbuf[1], qrsbuf, MEMMOVELEN);
                            qrsbuf[0] = newPeak;
                            qmedian = median(qrsbuf, 8);
                            det_thresh = thresh(qmedian, nmedian);
                            itaShifRight(rrbuf, MEMMOVELEN); //memmove(&rrbuf[1], rrbuf, MEMMOVELEN);
                            rrbuf[0] = count - QrsDetH.WINDOW_WIDTH;
                            rrmedian = median(rrbuf, 8);
                            sbcount = rrmedian + (rrmedian >> 1) + QrsDetH.WINDOW_WIDTH;
                            count = QrsDetH.WINDOW_WIDTH;

                            sbpeak = 0;

                            lastmax = maxder;
                            maxder = 0;
                            QrsDelay = QrsDetH.WINDOW_WIDTH + QrsDetH.FILTER_DELAY;
                            initBlank = initMax = rsetCount = 0;

                            // preBlankCnt = QrsDetH.PRE_BLANK ;
                        }
                        // If a peak isn't a QRS update noise buffer and estimate.
                        // Store the peak for possible search back.
                        else
                        {
                            itaShifRight(noise, MEMMOVELEN); //memmove(&noise[1], noise, MEMMOVELEN);
                            noise[0] = newPeak;
                            nmedian = median(noise, 8);
                            det_thresh = thresh(qmedian, nmedian);

                            // Don't include early peaks (which might be T-waves)
                            // in the search back process.  A T-wave can mask
                            // a small following QRS.

                            if ((newPeak > sbpeak) && ((count - QrsDetH.WINDOW_WIDTH) >= QrsDetH.MS360))
                            {
                                sbpeak = newPeak;
                                sbloc = count - QrsDetH.WINDOW_WIDTH;
                            }
                        }
                    }
                }
    		
                // Test for search back condition.  If a QRS is found in
                // search back update the QRS buffer and det_thresh.  

                if ((count > sbcount) && (sbpeak > (det_thresh >> 1)))
                {
                    itaShifRight(rrbuf, MEMMOVELEN); //memmove(&qrsbuf[1], qrsbuf, MEMMOVELEN);
                    qrsbuf[0] = sbpeak;
                    qmedian = median(qrsbuf, 8);
                    det_thresh = thresh(qmedian,nmedian);
                    itaShifRight(rrbuf, MEMMOVELEN); //memmove(&rrbuf[1], rrbuf, MEMMOVELEN);
                    rrbuf[0] = sbloc;
                    rrmedian = median(rrbuf, 8);
                    sbcount = rrmedian + (rrmedian >> 1) + QrsDetH.WINDOW_WIDTH;
                    QrsDelay = count = count - sbloc;
                    QrsDelay += QrsDetH.FILTER_DELAY;
                    sbpeak = 0;
                    lastmax = maxder;
                    maxder = 0;
                    initBlank = initMax = rsetCount = 0;
                }
            }

            // In the background estimate threshold to replace adaptive threshold
            // if eight seconds elapses without a QRS detection.

            if (qpkcnt == 8)
            {
                if (++initBlank == QrsDetH.MS1000)
                {
                    initBlank = 0;
                    rsetBuff[rsetCount] = initMax;
                    initMax = 0;
                    ++rsetCount;

                    // Reset threshold if it has been 8 seconds without
                    // a detection.

                    if(rsetCount == 8)
                    {
                        for(i = 0; i < 8; ++i)
                        {
                            qrsbuf[i] = rsetBuff[i];
                            noise[i] = 0;
                        }
                        qmedian = median(rsetBuff, 8);
                        nmedian = 0;
                        rrmedian = QrsDetH.MS1000;
                        sbcount = QrsDetH.MS1500 + QrsDetH.MS150;
                        det_thresh = thresh(qmedian, nmedian);
                        initBlank = initMax = rsetCount = 0;
                        sbpeak = 0;
                    }
                }
                if (newPeak > initMax)
                    initMax = newPeak;
            }

            return (QrsDelay);
	    }


        /***************************************************************************
        * peak() takes a datum as input and returns a peak height
        * when the signal returns to half its peak height, or 
        ****************************************************************************/

        unsafe public static int Peak(int datum, int init)
	    {
            int max = 0, timeSinceMax = 0, lastDatum = 0;
            int pk = 0;

            if (init != 0)
                max = timeSinceMax = 0;

            if (timeSinceMax > 0)
                ++timeSinceMax;

            if ((datum > lastDatum) && (datum > max))
            {
                max = datum;
                if (max > 2)
                    timeSinceMax = 1;
            }

            else if (datum < (max >> 1))
            {
                pk = max;
                max = 0;
                timeSinceMax = 0;
                Dly = 0;
            }

            else if (timeSinceMax > QrsDetH.MS95)
            {
                pk = max;
                max = 0;
                timeSinceMax = 0;
                Dly = 3;
            }
            lastDatum = datum;
            return (pk);
	    }


        /***************************************************************************
        median returns the median of an array of integers.  It uses a slow
        sort algorithm, but these arrays are small, so it hardly matters.
        ****************************************************************************/

        unsafe public static int median(int[] array, int datnum)
	    {
            int i, j, k, temp;
            int[] sort = new int[20];
            for(i = 0; i < datnum; ++i)
                sort[i] = array[i];

            for(i = 0; i < datnum; ++i)
            {
                temp = sort[i];
                for (j = 0; (temp < sort[j]) && (j < i); ++j) { }
                for (k = i - 1; k >= j; --k)
                    sort[k+1] = sort[k];
                sort[j] = temp;
            }
            return (sort[datnum >> 1]);

	    
        /*
        int median(int *array, int datnum)
        {
        long sum ;
        int i ;

        for(i = 0, sum = 0; i < datnum; ++i)
	        sum += array[i] ;
        sum /= datnum ;
        return(sum) ;
        } */
        } 


        /****************************************************************************
         thresh() calculates the detection threshold from the qrs median and noise
         median estimates.
        ****************************************************************************/

        unsafe public static int thresh(int qmedian, int nmedian)
	    {
            int thrsh, dmed;
            double temp;
            dmed = qmedian - nmedian;
            //	thrsh = nmedian + (dmed>>2) + (dmed>>3) + (dmed>>4); 
            temp = dmed;
            temp *= TH;
            dmed = (int) temp;
            thrsh = nmedian + dmed; // dmed * THRESHOLD 
            return (thrsh);
	    }


        /****************************************************************************
	        BLSCheck() reviews data to see if a baseline shift has occurred.
	        This is done by looking for both positive and negative slopes of
	        roughly the same magnitude in a 220 ms window.
        *****************************************************************************/

        unsafe public static int BLSCheck(int[] dBuf, int dbPtr, ref int maxder)
	    {
            int max, min;
            int maxt = 0, mint = 0;
            int t, x;
            max = min = 0;

            //return (0);

            for (t = 0; t < QrsDetH.MS220; ++t)
            {
                x = dBuf[dbPtr];
                if (x > max)
                {
                    maxt = t;
                    max = x;
                }
                else if (x < min)
                {
                    mint = t;
                    min = x;
                }
                if (++dbPtr == QrsDetH.DER_DELAY)
                    dbPtr = 0;
            }

            maxder = max;
            min = -min;

            // Possible beat if a maximum and minimum pair are found where the interval between them is less than 150 ms.

            if ((max > (min >> 3)) && (min > (max >> 3)) && (Math.Abs(maxt - mint) < QrsDetH.MS150))
                return (0);
            else
                return (1);
	    }


        //remake
        public static void itaShifRight(int[] list, int num)
        {
            //memmove(&rrbuf[1], rrbuf, MEMMOVELEN);
            num /= sizeof(int);
            for (int i = num; i > 0; i--)
                list[i] = list[i - 1];
        }
    }


}
